Rostagno Andrea
295706
Stima del reddito annuale in base ai dati del censimento¶
Introduzione¶
In questo elaborato presenteremo ed analizzeremo il dataset "Adult" riguardante il reddito medio annuale di un campione di persone in base ai dati del censimento.
I dati sono stati presi seguendo alcune regole come:
- l'età deve essere superiore ai 16 anni
- il reddito lordo rettificato della persona deve essere superiore a 100
- il "peso" di ogni dato deve essere maggiore di 1
- La persona deve lavorare più di 0 ore alla settimana
- Il compito di previsione consiste nel determinare se una persona guadagna oltre 50.000$ all'anno.
L'analisi procederà in questa maniera:
- Pulizia del database e visualizzazione di alcuni dati
- Utilizzo di tecniche per la riduzione di dimensionalità messe a confronto come la PCA e la MDA
- Utilizzo di tecniche di predizione come la LDA e la SVM
Durante il procedimento cercheremo di visualizzare il più possibile i vari risultati ottenuti.
Il dataset¶
Qui in allegato è presente il link diretto del dataset:
https://archive.ics.uci.edu/dataset/2/adult
Il dataset contiene 32561 campioni relativi a persone differenti e 14 features che andremo ora ad elencare.
- Age (Età): assume valori interi e non presenta valori mancanti
- Workclass (Classe di lavoro): Private, Self-emp-not-inc, Self-emp-inc, Federal-gov, Local-gov, State-gov, Without-pay, Never-worked. Presenta valori mancanti
- Fnlwgt (Peso finale): indica quella data persona quanta popolazione rappresenta, serve per fare una media ponderata. Assume solo valori interi e non presenta valori mancanti
- Education (Titoli di studio): Bachelors, Some-college, 11th, HS-grad, Prof-school, Assoc-acdm, Assoc-voc, 9th, 7th-8th, 12th, Masters, 1st-4th, 10th, Doctorate, 5th-6th, Preschool. Non presenta dati mancanti
- Education-num (Voto): Solo numeri interi e nessun dato mancante
- Marital-status (Stato civile): Married-civ-spouse, Divorced, Never-married, Separated, Widowed, Married-spouse-absent, Married-AF-spouse. Nessun dato mancante
- Occupation (Lavoro): Tech-support, Craft-repair, Other-service, Sales, Exec-managerial, Prof-specialty, Handlers-cleaners, Machine-op-inspct, Adm-clerical, Farming-fishing, Transport-moving, Priv-house-serv, Protective-serv, Armed-Forces. Sono presenti dei dati mancanti
- Relationship (Relazione): Wife, Own-child, Husband, Not-in-family, Other-relative, Unmarried. Nessun dato mancante
- Race(Etnia): White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other, Black. Nessun dato mancante
- Sex(Sesso): Female, Male. Nessun dato mancante
- Capital-gain(): Rappresenta il guadagno di capitale, ovvero quando vendi un qualcosa ad un prezzo superiore del prezzo di acquisto (p.e. acquisti delle azioni a 1000$ e le rivendi a 1500$, il guadagno capitale è di 500$). Assume solo valori interi e non presenta nessun dato mancante.
- Capital-loss(): Rappresenta la perdita di capitale, ovvero quando vendi un qualcosa ad un prezzo inferiore del prezzo di acquisto (p.e. acquisti delle azioni a 1500$ e le rivendi a 1000$, la perdita capitale è di 500$). Assume solo valori interi e non presenta nessun dato mancante.
- Hours-per-week(Ore settimanali lavorative): Assume solo valori interi e non presenta nessun dato mancante.
- Native-country(Nazionalità): United-States, Cambodia, England, Puerto-Rico, Canada, Germany, Outlying-US(Guam-USVI-etc), India, Japan, Greece, South, China, Cuba, Iran, Honduras, Philippines, Italy, Poland, Jamaica, Vietnam, Mexico, Portugal, Ireland, France, Dominican-Republic, Laos, Ecuador, Taiwan, Haiti, Columbia, Hungary, Guatemala, Nicaragua, Scotland, Thailand, Yugoslavia, El-Salvador, Trinadad&Tobago, Peru, Hong, Holand-Netherlands. Presenta dei dati mancanti
Invece per quanto riguarda gli attributi avremo la casella income (guadagno) che avrà come valori >50K, <=50K.
1. Pulizia e visualizzazione del dataset¶
%matplotlib widget
import pandas as pd
import numpy as np
from sklearn import datasets
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from FisherDA import MultipleFisherDiscriminantAnalysis as MDA
from sklearn.decomposition import PCA
from sklearn.svm import LinearSVC
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import display
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import precision_score, recall_score, fbeta_score
#importiamo il dataset
adult = pd.read_csv('adult.data', header=None,na_values=" ?")
adult_columns = ['age', 'workclass', 'fnlwgt', 'education', 'education-num', 'marital-status', 'occupation', 'relationship','race','sex','capital-gain','capital-loss','hours-per-week','native-country','income']
adult.rename(columns={k:adult_columns[k] for k in range(len(adult_columns))}, inplace=True)
display(adult)
| age | workclass | fnlwgt | education | education-num | marital-status | occupation | relationship | race | sex | capital-gain | capital-loss | hours-per-week | native-country | income | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 39 | State-gov | 77516 | Bachelors | 13 | Never-married | Adm-clerical | Not-in-family | White | Male | 2174 | 0 | 40 | United-States | <=50K |
| 1 | 50 | Self-emp-not-inc | 83311 | Bachelors | 13 | Married-civ-spouse | Exec-managerial | Husband | White | Male | 0 | 0 | 13 | United-States | <=50K |
| 2 | 38 | Private | 215646 | HS-grad | 9 | Divorced | Handlers-cleaners | Not-in-family | White | Male | 0 | 0 | 40 | United-States | <=50K |
| 3 | 53 | Private | 234721 | 11th | 7 | Married-civ-spouse | Handlers-cleaners | Husband | Black | Male | 0 | 0 | 40 | United-States | <=50K |
| 4 | 28 | Private | 338409 | Bachelors | 13 | Married-civ-spouse | Prof-specialty | Wife | Black | Female | 0 | 0 | 40 | Cuba | <=50K |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 32556 | 27 | Private | 257302 | Assoc-acdm | 12 | Married-civ-spouse | Tech-support | Wife | White | Female | 0 | 0 | 38 | United-States | <=50K |
| 32557 | 40 | Private | 154374 | HS-grad | 9 | Married-civ-spouse | Machine-op-inspct | Husband | White | Male | 0 | 0 | 40 | United-States | >50K |
| 32558 | 58 | Private | 151910 | HS-grad | 9 | Widowed | Adm-clerical | Unmarried | White | Female | 0 | 0 | 40 | United-States | <=50K |
| 32559 | 22 | Private | 201490 | HS-grad | 9 | Never-married | Adm-clerical | Own-child | White | Male | 0 | 0 | 20 | United-States | <=50K |
| 32560 | 52 | Self-emp-inc | 287927 | HS-grad | 9 | Married-civ-spouse | Exec-managerial | Wife | White | Female | 15024 | 0 | 40 | United-States | >50K |
32561 rows × 15 columns
adult.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 32561 entries, 0 to 32560 Data columns (total 15 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 age 32561 non-null int64 1 workclass 30725 non-null object 2 fnlwgt 32561 non-null int64 3 education 32561 non-null object 4 education-num 32561 non-null int64 5 marital-status 32561 non-null object 6 occupation 30718 non-null object 7 relationship 32561 non-null object 8 race 32561 non-null object 9 sex 32561 non-null object 10 capital-gain 32561 non-null int64 11 capital-loss 32561 non-null int64 12 hours-per-week 32561 non-null int64 13 native-country 31978 non-null object 14 income 32561 non-null object dtypes: int64(6), object(9) memory usage: 3.7+ MB
#elimino tutte le righe con valori nulli e salvo le dimensioni
adult_b = adult.dropna()
M,N = adult_b.shape
adult['income'].value_counts()
income <=50K 24720 >50K 7841 Name: count, dtype: int64
#mostro che la maggior parte dei dati eliminati guadagna meno di 50k(però non è un problema perchè la differenza è ancora sostanziale)
adult_b['income'].value_counts()
income <=50K 22654 >50K 7508 Name: count, dtype: int64
#una distribuzione dell'età
plt.figure(figsize=(8, 4))
plt.hist(adult_b['age'], bins=25, edgecolor='black')
plt.title('Distribuzione dell\'età')
plt.xlabel('Età')
plt.ylabel('Frequenza')
plt.xticks(ticks=np.arange(16,92,5), labels=[f'{i}' for i in range(16,92,5)])
plt.yticks(ticks=np.arange(0,2750,250), labels=[f'{i}' for i in range(0,2750,250)])
plt.grid(False)
plt.show()
average_age = np.round((adult_b['age']*adult_b['fnlwgt']).sum()/adult_b['fnlwgt'].sum())
display("L'età media ponderata è: "+str(average_age)+" anni")
"L'età media ponderata è: 38.0 anni"
#trasformo le colonne object in category
adult_b=adult_b.astype({"workclass":'category', "education":'category', "marital-status":'category', "occupation":'category', "relationship":'category',"race":'category',"sex":'category',"native-country":'category',"income":'category'})
#Visto che le classi income sono due (>50k, <=50k), per esse possiamo anche utilizzare un label encoding. Per le features adotteremo invece sempre il one-hot encoding
adult_work=adult_b.copy()
le=LabelEncoder()
le.fit(adult_work['income'].values)
adult_work['income']=le.transform(adult_work['income'].values)
#one-hot encoding
onehot_cols=['workclass', 'education', 'marital-status', 'occupation', 'relationship', 'race','sex','native-country']
adult_work_old=adult_work.copy()
adult_work=pd.get_dummies(adult_work, columns=onehot_cols)
adult_work=adult_work.astype(int)
col_to_move = adult_work.columns[6]
new_pos = 104
cols = adult_work.columns.tolist()
cols.insert(new_pos, cols.pop(cols.index(col_to_move)))
adult_work = adult_work[cols]
display(adult_work)
| age | fnlwgt | education-num | capital-gain | capital-loss | hours-per-week | workclass_ Federal-gov | workclass_ Local-gov | workclass_ Private | workclass_ Self-emp-inc | ... | native-country_ Puerto-Rico | native-country_ Scotland | native-country_ South | native-country_ Taiwan | native-country_ Thailand | native-country_ Trinadad&Tobago | native-country_ United-States | native-country_ Vietnam | native-country_ Yugoslavia | income | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 39 | 77516 | 13 | 2174 | 0 | 40 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 1 | 50 | 83311 | 13 | 0 | 0 | 13 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 | 38 | 215646 | 9 | 0 | 0 | 40 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 3 | 53 | 234721 | 7 | 0 | 0 | 40 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 4 | 28 | 338409 | 13 | 0 | 0 | 40 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 32556 | 27 | 257302 | 12 | 0 | 0 | 38 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 32557 | 40 | 154374 | 9 | 0 | 0 | 40 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 |
| 32558 | 58 | 151910 | 9 | 0 | 0 | 40 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 32559 | 22 | 201490 | 9 | 0 | 0 | 20 | 0 | 0 | 1 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| 32560 | 52 | 287927 | 9 | 15024 | 0 | 40 | 0 | 0 | 0 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 |
30162 rows × 105 columns
#mostriamo la relazione tra età e ore settimanali lavorative
cmap=ListedColormap(['blue', 'orange'])
plt.figure(figsize=(8, 4))
plt.scatter(adult_work_old['age'], adult_work_old['hours-per-week'], c=adult_work_old['income'],cmap=cmap, alpha=0.5)
plt.title('Relazione tra Età e ore lavorative settimanali')
plt.xlabel('Età')
plt.ylabel('ore lavorative settimanali')
handles = [plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=cmap(0), markersize=10, label='<=50K'),
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=cmap(1), markersize=10, label='>50K')]
plt.legend(handles=handles, title='income')
plt.xticks(ticks=np.arange(16,92,5), labels=[f'{i}' for i in range(16,92,5)])
plt.yticks(ticks=np.arange(0,110,10), labels=[f'{i}' for i in range(0,110,10)])
plt.grid(True)
plt.show()
Notiamo che tra i 16 e i 26 anni tutte le persone che lavorano fino a 40 ore settimanali non guadagnano più di 50k al mese.
Le persone che guadagnano di più invece si trovano in un età compresa tra i 31 e i 61 anni e lavorano tra le 40 e le 70 ore a settimana.
#mostriamo la relazione tra guadagni da investimenti ed età
plt.figure(figsize=(8, 4))
plt.scatter(adult_work_old['age'], adult_work_old['capital-gain'], c=adult_work_old['income'], cmap=cmap, alpha=0.50)
plt.title('Età vs Guadagno di Capitale')
plt.xlabel('Età')
plt.ylabel('Guadagno di Capitale')
plt.legend(handles=handles, title='Income')
plt.xticks(ticks=np.arange(16,92,5), labels=[f'{i}' for i in range(16,92,5)])
plt.yticks(ticks=np.arange(0,110000,10000), labels=[f'{i}' for i in range(0,110000,10000)])
plt.grid(True)
plt.ion()
plt.show()
Anche qui un altro dato interessante, praticamente tutti quelli che guadagnano più di 50k all'anno percepiscono almeno 10k da investimenti.
Mentre chi guadagna meno di 50k non percepisce oltre 5.000$ da investimenti vari.
#mostriamo la relazione tra ore lavorative settimanali e il peso che hanno sui dati
plt.figure(figsize=(8, 4))
plt.scatter(adult_work_old['age'], adult_work_old['capital-loss'], c=adult_work_old['income'], cmap=cmap, alpha=0.5)
plt.title('Età vs Perdita di Capitale')
plt.xlabel('Età')
plt.ylabel('Perdita di Capitale')
plt.legend(handles=handles, title='Income')
plt.xticks(ticks=np.arange(16,92,5), labels=[f'{i}' for i in range(16,92,5)])
plt.yticks(ticks=np.arange(0,5000,500), labels=[f'{i}' for i in range(0,5000,500)])
plt.grid(True)
plt.show()
Chi guadagna di più sembra avere due "stop-loss", il primo a 2000 e il secondo a 2500.
Si nota inoltre che le perdite per investimenti sono concentrate tra i 16 e i 61 anni di età, successivamente diminuiscono nettamente.
#mostriamo la correlazione tra le variabili
adult_work.iloc[:, :-1].corr()
| age | fnlwgt | education-num | capital-gain | capital-loss | hours-per-week | workclass_ Federal-gov | workclass_ Local-gov | workclass_ Private | workclass_ Self-emp-inc | ... | native-country_ Portugal | native-country_ Puerto-Rico | native-country_ Scotland | native-country_ South | native-country_ Taiwan | native-country_ Thailand | native-country_ Trinadad&Tobago | native-country_ United-States | native-country_ Vietnam | native-country_ Yugoslavia | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age | 1.000000 | -0.076511 | 0.043526 | 0.080154 | 0.060165 | 0.101599 | 0.056626 | 0.068256 | -0.210491 | 0.111039 | ... | 0.001062 | 0.007836 | 0.000024 | 0.001923 | -0.007879 | -0.004940 | 0.007868 | 0.016259 | -0.017775 | 0.000657 |
| fnlwgt | -0.076511 | 1.000000 | -0.044992 | 0.000422 | -0.009750 | -0.022886 | -0.006932 | -0.003520 | 0.046589 | -0.025496 | ... | -0.014020 | 0.007121 | -0.003055 | -0.010598 | 0.001615 | -0.001241 | -0.000413 | -0.083390 | -0.010761 | 0.005707 |
| education-num | 0.043526 | -0.044992 | 1.000000 | 0.124416 | 0.079646 | 0.152522 | 0.058244 | 0.097378 | -0.165069 | 0.078843 | ... | -0.043058 | -0.042083 | 0.001815 | 0.017812 | 0.049129 | 0.008183 | -0.017134 | 0.127207 | -0.010953 | -0.001661 |
| capital-gain | 0.080154 | 0.000422 | 0.124416 | 1.000000 | -0.032229 | 0.080432 | -0.006299 | -0.009624 | -0.048185 | 0.096482 | ... | -0.003826 | -0.006270 | -0.002816 | -0.002582 | 0.007639 | -0.003501 | -0.003603 | 0.012375 | -0.002493 | -0.002317 |
| capital-loss | 0.060165 | -0.009750 | 0.079646 | -0.032229 | 1.000000 | 0.052417 | 0.010380 | 0.014727 | -0.036377 | 0.030956 | ... | -0.007343 | -0.004560 | -0.004175 | 0.005677 | 0.005679 | -0.005191 | 0.008849 | 0.015119 | 0.000344 | -0.005036 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| native-country_ Thailand | -0.004940 | -0.001241 | 0.008183 | -0.003501 | -0.005191 | 0.012846 | -0.004266 | -0.006441 | -0.008144 | 0.018052 | ... | -0.000798 | -0.001430 | -0.000454 | -0.001154 | -0.000887 | 1.000000 | -0.000580 | -0.076390 | -0.001095 | -0.000547 |
| native-country_ Trinadad&Tobago | 0.007868 | -0.000413 | -0.017134 | -0.003603 | 0.008849 | -0.007566 | -0.004390 | -0.001255 | 0.002164 | 0.002630 | ... | -0.000821 | -0.001472 | -0.000467 | -0.001187 | -0.000913 | -0.000580 | 1.000000 | -0.078606 | -0.001127 | -0.000563 |
| native-country_ United-States | 0.016259 | -0.083390 | 0.127207 | 0.012375 | 0.015119 | 0.010673 | 0.017541 | 0.032937 | -0.049809 | 0.007350 | ... | -0.108062 | -0.193727 | -0.061442 | -0.156254 | -0.120121 | -0.076390 | -0.078606 | 1.000000 | -0.148334 | -0.074108 |
| native-country_ Vietnam | -0.017775 | -0.010761 | -0.010953 | -0.002493 | 0.000344 | -0.010381 | -0.000004 | -0.003952 | 0.009369 | -0.004973 | ... | -0.001549 | -0.002777 | -0.000881 | -0.002240 | -0.001722 | -0.001095 | -0.001127 | -0.148334 | 1.000000 | -0.001062 |
| native-country_ Yugoslavia | 0.000657 | 0.005707 | -0.001661 | -0.002317 | -0.005036 | 0.006983 | -0.004139 | -0.000550 | 0.003861 | 0.003343 | ... | -0.000774 | -0.001387 | -0.000440 | -0.001119 | -0.000860 | -0.000547 | -0.000563 | -0.074108 | -0.001062 | 1.000000 |
104 rows × 104 columns
PCA¶
Introduzione
Quando ci troviamo ad operare in situazioni in cui si hanno dati con un gran numero di dimensioni, è naturale pensare di
proiettare questi dati in un sottospazio di dimensionalità inferiore, cercando di non perdere informazione importante sulle variabili originarie.
Un modo per ottenere tale riduzione è attraverso la selezione delle variabili, chiamata anche “feature selection”.
La Principal Component Analysis (PCA) è una tecnica finalizzata a ridurre la dimensionalità di un insieme di dati con
finalità esplorativa o di visualizzazione dei dati, per un eventuale uso in analisi successive.
Procedimento
Supponiamo di trovarci in uno spazio $\mathbb{R}^{d}$ e di volere una funzione lineare ${f}$ tali che ${f} : \mathbb{R}^d\rightarrow \mathbb{R}^e$ ma al contempo un'altra funzione lineare ${g}$ tali che ${g} : \mathbb{R}^e\rightarrow\mathbb{R}^d$, con ${e}$ < ${d}$ .
Associamo a queste due trasformazioni le loro corrispettive matrici ${U}\in\mathbb{R}^{e\times d}$ e ${V}\in\mathbb{R}^{d\times e}$
Il nostro campione ${S}$ ha distribuzione ${D}^{m}$ ovvero ${S}=\{\boldsymbol{p}_1,\ldots ,\boldsymbol{p}_m\}$.
Il nostro obbiettivo è quello di trovare le due migliori funzioni lineari che permettano di fare questo procedimento minimizzando la perdita di informazione.
Formalizziamo:
$\underset{U,V}{\text{argmin}} \sum_{i=1}^{m} \left\| p_i - {VU}p_i \right\|^2_2$
Ora dobbiamo trovare una matrice ${U}\in\mathbb{R}^{e\times d}$ e una matrice ${V}\in\mathbb{R}^{d\times e}$ che siano soluzione del problema di minimizzazione.
- Dobbiamo calcolare una matrice ${A}=\sum_{i=1}^{m} \vec{p_i} \vec{p_i}^T$ dove $p_i$ è un vettore noto di dimensione ${d} \times{1}$.
- La matrice ${A}$ può anche essere calcolata come ${A}={X}^T{X}$ dove ${X}$ è la matrice che ha come righe i campioni di ${S}$ quindi avrà dimensione ${m} \times{d}$.
- Notiamo inoltre che ${A}^T={X}^T{X}={A}$ quindi la matrice ${A}$ è semidefinita positiva.
- Notiamo inoltre che ${A}^T={X}^T{X}={A}$ quindi la matrice ${A}$ è semidefinita positiva.
- Tramite la scomposizione (SVD) diventa ${A}={B}{D}{B}^T$ dove ${D}$ è la matrice con gli autovalori di ${A}$ posti in ordine decrescente mentre la matrice ${B}$ è la matrice con colonne gli autovettori ad essi associati.
- Possiamo dedurre che le componenti principali saranno gli autovettori della matrice ${B}$ e che le matrici incognite saranno appunto ${V}$ con colonne i primi ${e}$ autovettori e ${U}$ sarà uguale a ${V}^T$.
- Possiamo dedurre che le componenti principali saranno gli autovettori della matrice ${B}$ e che le matrici incognite saranno appunto ${V}$ con colonne i primi ${e}$ autovettori e ${U}$ sarà uguale a ${V}^T$.
Riformuliamo il problema:
$\underset{V}{\text{argmax}} \text{ trace} ({V}^T(\sum_{i=1}^m \vec{p_i} \vec{p_i}^T ){V}) $
${V}{V}^T={I}$
${V}\in\mathbb{R}^{d\times e}$
Adesso che abbiamo ottenuto le matrice ${V}$ (${U}={V}^T$ ) e che sappiamo che ha come colonne gli autovettori associati agli autovalori di ${A}$ in ordine decrescente ci domandiamo:
ma quanto stiamo sbagliando a passre da uno spazio di dimensione ${d}$ ad uno spazio di dimensione ${e}$ ?
$\sum_{i={e}+1}^{d} \lambda_i({A}) $
Ovvero sommiamo gli autovalori che sono stati eliminati.
Varianza spiegata
Per costruzione, la varianza totale si può scrivere: $\sum_{i=1}^{d} \lambda_i = \text{tr}({A})$
dove il parametro $\lambda_i$ rappresenta la varianza spiegata dall' i-esima componente principale.
$\frac{\lambda_i}{\text{tr}({A})} = \frac{\lambda_i}{\sum_{j}\lambda_j} = \text{PVE}_i$
Grazie a questo calcolo si può quantificare quanta "informazione" porta ogni componente principale.
Non c'è un modo preciso per scegliere quante componenti principali vogliamo, ma generalmente si sceglie in base a quanta informazione vogliamo tenere (p.e. 0.90 , quante PC corrispondono? ) o a problemi di rappresentazione (p.e. 2 o 3 PC per la raffigurazione dei dati).
Un altro modo di visualizzare la varianza è tramite la varianza cumulata (CPVE):
a differenza della PVE che calcola in ogni PC la porzione di varianza che contiene, la CPVE mostra quanta varianza è stata rappresentata fino a quella PC.
Sia ${k}$ riferito alla PC di riferimento, allora la CPVE in quel punto sarà: $\frac{\sum_{i=1}^{k}\lambda_i}{\text{tr}({A})} = \text{CPVE}_i$
#separiamo le classi dalle features
x=adult_work.iloc[:, :-1].values
y=adult_work['income'].values
#generiamo il training
random_seed = 42
test_p = 0.50
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=test_p, random_state=random_seed, shuffle=True)
#standardizzo
scaler_X = StandardScaler()
scaler_X.fit(X_train)
X_scaled = scaler_X.transform(X_train)
#applichiamo la PCA sia st che nonst
pca=PCA()
pca_nonst=PCA()
pca.fit(X_scaled)
pca_nonst.fit(X_train)
#mostriamo la varianza spiegata
plt.figure()
plt.plot(np.insert(np.cumsum(pca.explained_variance_ratio_), 0, 0))
plt.title('x_scaled')
plt.xticks(ticks=np.arange(1,pca.n_components_+3,8),
labels=[f'PC{i}' for i in range(1, pca.n_components_ + 3,8)], rotation=45)
plt.xlabel('Principal components')
plt.ylabel('Cumulative explained variance (%)')
plt.grid()
plt.show()
plt.figure()
plt.bar(np.arange(pca_nonst.n_components_), pca_nonst.explained_variance_ratio_, bottom=np.insert(np.cumsum(pca_nonst.explained_variance_ratio_), 0, 0)[:-1])
plt.plot(np.cumsum(pca_nonst.explained_variance_ratio_), 'r')
plt.title('x_nonst')
plt.xticks(ticks=np.arange(1,pca_nonst.n_components_,10),
labels=[f'PC{i}' for i in range(1, pca_nonst.n_components_ + 1,10)])
plt.xlabel('Principal components')
plt.ylabel('Cumulative explained variance (%)')
plt.grid()
plt.show()
#scegliamo il numero di componenti principali in base alla varianza spiegata
explvar_p = 0.90
pca=PCA(explvar_p)
pca.fit(X_scaled)
pca_adult=pca.transform(X_scaled)
#trasformo i miei dati originali in quelli trasformati dalla PCA
features_scalate=scaler_X.transform(x)
Zpca = pca.transform(features_scalate)
df_pca = pd.DataFrame({'val': [pca.n_components_, np.round(pca.explained_variance_ratio_.sum()*100,decimals=2)]}, index=['n. PC', 'expl. Var. (%)'])
display(df_pca)
| val | |
|---|---|
| n. PC | 79.00 |
| expl. Var. (%) | 90.51 |
#mostriamo il grafico a dispersione
plt.figure(figsize=(12, 12))
plt.scatter(pca_adult[:, 0], pca_adult[:, 1], c=y_train, cmap=cmap, alpha=0.5)
plt.title('ADULT - SCORE GRAPH - 2D')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(handles=handles, title='income')
plt.grid()
plt.show()
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(pca_adult[:, 0], pca_adult[:, 1],pca_adult[:,2], c=y_train, cmap=cmap, alpha=0.5)
plt.title('ADULT - SCORE GRAPH - 3D')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
plt.grid()
plt.show()
#mostriamo la correlazione tra le componenti principali e le features
plt.figure(figsize=(18, 15))
for i in range(pca.n_components_):
plt.plot([0, pca.components_[0, i]], [0, pca.components_[1, i]], label=adult_work.iloc[:, :-1].columns[i])
plt.scatter(pca.components_[0, i], pca.components_[1,i], c='k')
plt.legend(bbox_to_anchor=(1.05, 1), fontsize='x-small')
plt.title('ADULT - LOADING GRAPH 2D')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.grid()
plt.show()
fig = plt.figure(figsize=(18, 15))
ax = fig.add_subplot(111, projection='3d')
for i in range(pca.n_components_):
ax.plot([0, pca.components_[0, i]], [0, pca.components_[1, i]], [0, pca.components_[2, i]])
ax.scatter(pca.components_[0, i], pca.components_[1, i], pca.components_[2, i], c='k')
plt.title('ADULT - LOADING GRAPH 3D')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
plt.grid()
plt.show()
#tramite dei barplot mostriamo la correlazione tra le componenti principali e le features
plt.figure(figsize=(18, 15))
plt.bar(np.arange(0,pca.n_components_,1), pca.components_[0, :-25])
col_n=[None]*pca.n_components_
for i in range(pca.n_components_):
col_n[i]=adult_work.iloc[:, :-1].columns[i]
plt.xticks(ticks=np.arange(pca.n_components_), labels=col_n,rotation=90, fontsize='x-small')
plt.title('PC1')
plt.grid()
plt.show()
plt.figure(figsize=(18, 15))
plt.bar(np.arange(0,pca.n_components_,1), pca.components_[1, :-25])
col_n=[None]*pca.n_components_
for i in range(pca.n_components_):
col_n[i]=adult_work.iloc[:, :-1].columns[i]
plt.xticks(ticks=np.arange(pca.n_components_), labels=col_n,rotation=90, fontsize='x-small')
plt.title('PC2')
plt.grid()
plt.show()
MDA¶
Introduzione
Fisher Discriminant Analysis (FDA) è una tecnica di riduzione della dimensionalità e classificazione che proietta i dati suddivisi in $c\in\mathbb{N}$ classi, in uno spazio di dimensione $1\leq m \leq (c-1)$ cercando di massimizzare la separazione tra le classi dei dati.
L'obiettivo è trovare un asse che massimizzi la distanza tra i centri delle classi e minimizzi la varianza all'interno delle classi.
Usare la FDA equivale a dire che si usa una trasformazione che preserva la classificazione dei dati.
MDA è particolarmente utile quando si lavora con dataset di alta dimensione e si desidera ridurre la dimensionalità mantenendo la separazione tra le classi per migliorare la classificazione
(p.e. se nel dataset sono presenti solo due classi non è molto utile in quanto, rispetto alla formula, proietterà i dati in un'unica dimensione).
Differenze tra PCA e MDA¶
Principal Component Analysis (PCA) e Multiple Discriminant Analysis (MDA) (nota anche come Fisher Discriminant Analysis, FDA) sono entrambe tecniche di riduzione della dimensionalità, ma con obiettivi e approcci differenti.
PCA¶
- Obiettivo: Ridurre la dimensionalità massimizzando la varianza dei dati.
- Metodo: Trova le componenti principali che spiegano la maggior parte della varianza nei dati.
- Uso: Adatto per l'analisi esplorativa dei dati e la compressione dei dati.
- Unsupervised learning: Non utilizza informazioni sulle etichette delle classi.
MDA (FDA)¶
- Obiettivo: Ridurre la dimensionalità massimizzando la separazione tra le classi.
- Metodo: Trova la proiezione lineare che massimizza la distanza tra i centri delle classi e minimizza la varianza all'interno delle classi.
- Uso: Adatto per problemi di classificazione supervisionata.
- Supervised learning: Utilizza informazioni sulle etichette delle classi.
Procedimento
Utilizziamo un campione ${S}=\{\boldsymbol{x}_1,\ldots ,\boldsymbol{x}_k\}\subset\mathbb{R}^d$ e una vettore classe ${c}=\{\boldsymbol{c}_1,\ldots ,\boldsymbol{c}_3\}$ (supponiamo di avere 3 classi)
Calcolo dei Vettori Medi
Calcolare il vettore medio per ciascuna classe.
$\mu_i = \frac{1}{N_i} \sum_{x_k \in c_i} x_k$
dove ${N_i}$ è il numero di campioni della classe ${c_i}$ mentre ${x_k}$ sono sono i campioni di classe ${c_i}$.
Al termine di questa operazione avremo ${i}$ vettori media distinti (p.e. se ho 3 classi 0,1,2 avrò tre vettori) dove ciascuno sarà ottenuto esclusivamente dai campioni appartenenti a quella classe.Matrice di Scatter Within-Class (Sw)
Calcolare la matrice di scatter all'interno delle classi.
$S_w = \sum_{i=1}^{c} \sum_{x_k \in c_i} (x_k - \mu_i)(x_k - \mu_i)^T$
dove ${c}$ è il numero di classi.
Questa matrice $S_w$ somma tra loro le matrici ottenute centrando tutti i campioni con le loro rispettive medie e moltiplicandoli tra loro,
assomiglia a una parte della formula per calcolare la matrice di covarianza: $S = \frac{1}{N-1} \sum_{i=1}^{N} (x_i - \mu)(x_i - \mu)^T$ .Matrice di Scatter Between-Class (Sb)
Calcolare la matrice di scatter tra le classi.
$S_b = \sum_{i=1}^{c} N_i (\mu_i - \mu)(\mu_i - \mu)^T$
dove $\mu$ è il vettore medio globale ottenuto facendo $\frac{1}{k}\sum_{i=1}^{k}x_i$.
Autovettori e Autovalori
Risolvere il problema degli autovalori per $S_w^{-1}S_b$ per trovare gli autovettori (direzioni) e gli autovalori (importanza delle direzioni).
$S_w^{-1} S_b \vec{v} = \lambda \vec{v}$
dove $\vec{v}$ sono gli autovettori e $\lambda$ sono gli autovalori della matrice $S=S_w^{-1}S_b$.
Selezione delle Componenti Principali
Ordinare gli autovettori in base ai rispettivi autovalori in ordine decrescente e selezionare i primi ${n}$ autovettori per formare la matrice di trasformazione ${W}$ (uguale alla PCA).
Proiezione dei Dati
Proiettare i dati originali nello spazio delle componenti principali selezionate.
$Z = XW$
dove $Z$ sono i dati proiettati, $X$ è la matrice dei dati originali e $W$ è la matrice dei primi $n$ autovettori.
#mostriamo la separazione tra classi dell'MDA
mda=MDA()
mda.fit(X_train, y_train)
Zmda = mda.transform(x)
zero=np.zeros(Zmda.shape[0])
fig1, axs1 = plt.subplots(1, 2, figsize=(14, 10))
axs1[0].scatter(Zmda[:, 0],zero, c=y, cmap=cmap, alpha=0.50)
handles = [plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=cmap(0), markersize=10, label='<=50K'),
plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=cmap(1), markersize=10, label='>50K')]
axs1[0].legend(handles=handles, title='income')
axs1[0].set_title('MDA')
axs1[0].grid()
axs1[0].set_xticks(np.arange(Zmda[:, 0].min()-1, Zmda[:, 0].max()+1))
axs1[1].scatter(Zpca[:, 0], Zpca[:, 1], c=y, cmap=cmap, alpha=0.50)
axs1[1].legend(handles=handles, title='income')
axs1[1].set_title('PCA')
axs1[1].grid()
#applichiamo la MDA al campione della PCA
mda_pca=MDA()
mda_pca.fit(Zpca, y)
Zmda_pca = mda_pca.transform(Zpca)
fig1, axs1 = plt.subplots(1, 2, figsize=(14, 5))
axs1[0].scatter(Zmda[:,0].flatten(), y, c=y, cmap=cmap, alpha=0.15)
axs1[0].set_title('MDA (senza PCA)')
axs1[0].grid()
axs1[0].set_xticks(np.arange(Zmda[:, 0].min(), Zmda[:, 0].max()))
axs1[1].scatter(Zmda_pca.flatten(), y, c=y, cmap=cmap, alpha=0.15)
axs1[1].set_title('MDA (con PCA)')
axs1[1].grid()
LDA¶
Introduzione
Linear Discriminant Analysis (LDA) è una tecnica di classificazione supervisionata che cerca di trovare una combinazione lineare delle caratteristiche che meglio separa le classi (analogo come la MDA).
La novità è che durante il processo di predizione, LDA calcola le probabilità a priori delle classi e le probabilità condizionate dei dati rispetto a ciascuna classe.
Utilizzando il teorema di Bayes, LDA stima la probabilità posteriore di ogni classe data una nuova osservazione e assegna la classe con la massima probabilità posteriore, ottimizzando così la separazione e la classificazione dei dati.
Procedimento
Iniziamo con il mostrare il teorema di Bayes: $ P(Y=k \, |\, X=\boldsymbol{x}) = \frac{P(X=\boldsymbol{x} \, |\, Y=k)\cdot P(Y=k)}{P(X=\boldsymbol{x})}$
Calcolo delle Probabilità A Priori:
Sia ${\omega}=\{\boldsymbol{\omega}_1,\ldots ,\boldsymbol{\omega}_n\}$ il vettore contenente le diverse classi.
Calcolare la probabilità a priori $P(\omega_i)$ di ciascuna classe.$P(\omega_i) = \frac{N_i}{N}$
dove $N_i$ è il numero di campioni nella classe $ \omega_i $ e $ N $ è il numero totale di campioni.
Calcolo delle Probabilità Condizionate:
Per ciascuna classe, calcolare la probabilità condizionata $ P(x | \omega_i) $.
Poiché LDA assume una distribuzione normale multivariata per ciascuna classe, si utilizza la funzione di densità di probabilità della distribuzione normale:$P(x | \omega_i) = \frac{1}{(2\pi)^{d/2} |\Sigma|^{1/2}} \exp\left(-\frac{1}{2} (x - \mu_i)^T \Sigma^{-1} (x - \mu_i)\right)$
dove $ \mu_i $ è il vettore medio della classe $ \omega_i $, $ \Sigma $ è la matrice di covarianza (che è la stessa per tutte le classi in LDA) e $ d $ è la dimensione dei dati ${x}\in\mathbb{R}^{d}$.
Calcolo della Probabilità Posteriore:
Utilizzare il teorema di Bayes per calcolare la probabilità posteriore $ P(\omega_i | x) $:
$P(\omega_i | x) \propto P(x | \omega_i) P(\omega_i)$
Normalizzare per assicurarsi che le probabilità posteriori sommino a 1
(p.e. ho due classi 0,1 la prob. che il mio $x$ sia di classe 1 è 0.84 mentre che sia di classe 0 è del 0.16 . 0.84+0.16=1 ).Decisione di Classificazione:
Assegnare la classe che massimizza la probabilità posteriore:
$\hat{\omega} = \underset{\omega_i}{\text{argmax}} \, P(\omega_i | x)$
ovvero si assegna la classe che ha probabilità maggiore di appartenere a quel campione.
#calcolo la precisione con cui il metodo LDA predice dei dati(no PCA)
lda =LDA()
lda.fit(X_train, y_train)
y_pred=lda.predict(X_test)
y_pred_proba=lda.predict_proba(X_test)
y_pred_df= pd.DataFrame({'Pred. Class': y_pred,
'P(Class 0) - %': np.round(y_pred_proba[:, 0] * 100, decimals=2),
'P(Class 1) - %': np.round(y_pred_proba[:, 1] * 100, decimals=2)})
scores_dict = {'Training Set': np.round(lda.score(X_train, y_train)*100, decimals=2), 'Test Set': np.round(lda.score(X_test, y_test) * 100, decimals=2)}
scores = pd.DataFrame(scores_dict, index=['Accuracy'])
display(scores)
display(y_pred_df)
| Training Set | Test Set | |
|---|---|---|
| Accuracy | 83.79 | 83.78 |
| Pred. Class | P(Class 0) - % | P(Class 1) - % | |
|---|---|---|---|
| 0 | 0 | 64.07 | 35.93 |
| 1 | 0 | 85.07 | 14.93 |
| 2 | 0 | 63.54 | 36.46 |
| 3 | 0 | 81.56 | 18.44 |
| 4 | 0 | 90.70 | 9.30 |
| ... | ... | ... | ... |
| 15076 | 0 | 69.66 | 30.34 |
| 15077 | 0 | 64.27 | 35.73 |
| 15078 | 1 | 40.30 | 59.70 |
| 15079 | 0 | 89.76 | 10.24 |
| 15080 | 0 | 99.20 | 0.80 |
15081 rows × 3 columns
#calcolo la precisione del metodo LDA con PCA
lda_pca =LDA()
lda_pca.fit(pca_adult, y_train)
y_pred=lda_pca.predict(pca.transform(scaler_X.transform(X_test)))
y_pred_proba=lda_pca.predict_proba(pca.transform(scaler_X.transform(X_test)))
y_pred_df= pd.DataFrame({'Pred. Class': y_pred,
'P(Class 0) - %': np.round(y_pred_proba[:, 0] * 100, decimals=2),
'P(Class 1) - %': np.round(y_pred_proba[:, 1] * 100, decimals=2)})
scores_dict = {'Training Set': np.round(lda_pca.score(pca_adult, y_train)*100, decimals=2), 'Test Set': np.round(lda_pca.score(pca.transform(scaler_X.transform(X_test)), y_test) * 100, decimals=2)}
scores = pd.DataFrame(scores_dict, index=['Accuracy'])
display(scores)
display(y_pred_df)
| Training Set | Test Set | |
|---|---|---|
| Accuracy | 83.52 | 83.51 |
| Pred. Class | P(Class 0) - % | P(Class 1) - % | |
|---|---|---|---|
| 0 | 0 | 51.79 | 48.21 |
| 1 | 0 | 78.36 | 21.64 |
| 2 | 0 | 68.79 | 31.21 |
| 3 | 0 | 78.25 | 21.75 |
| 4 | 0 | 83.97 | 16.03 |
| ... | ... | ... | ... |
| 15076 | 0 | 71.38 | 28.62 |
| 15077 | 1 | 41.91 | 58.09 |
| 15078 | 1 | 40.30 | 59.70 |
| 15079 | 0 | 88.85 | 11.15 |
| 15080 | 0 | 98.34 | 1.66 |
15081 rows × 3 columns
SVM Lineare¶
Introduzione
Il problema della SVM lineare consiste nel trovare un iperpiano in $\mathbb{R}^n$ che separi i punti dati in due classi, una positiva e una negativa.
Un iperpiano è definito come l'insieme di punti x che soddisfano l'equazione $w · x + b = 0$, dove ${w}\in\mathbb{R}^{n}$ è un vettore non nullo normale all'iperpiano e ${b}\in\mathbb{R}$ è uno scalare.
L'iperpiano divide lo spazio in due semispazi.
Procedimento¶
Caso separabile
Si assume che esista un iperpiano che separi perfettamente i punti di training in due popolazioni di punti etichettati positivamente e negativamente.
Tuttavia, ci sono infiniti iperpiani che potrebbero separare i dati.
L'obiettivo delle SVM è trovare l'iperpiano che massimizza il margine, ovvero la distanza tra l'iperpiano e il punto dati più vicino.
Questo iperpiano è chiamato iperpiano a margine massimo e la distanza tra l'iperpiano e i punti dati più vicini è chiamata margine geometrico.
I punti dati che si trovano sul margine sono chiamati vettori di supporto.
Indichiamo rispettivamente l'insieme dei vettori e l'insieme delle classi con ${X}=\{\boldsymbol{x}_1,\ldots ,\boldsymbol{x}_n\}$ e ${Y}=\{\boldsymbol{y}_1,\ldots ,\boldsymbol{y}_n\}$
La funzione obbiettivo diventa:
$
\begin{aligned}
& \underset{w, b}{\text{min}} \, \frac{1}{2} \|w\|^2 \\
& y_i (w \cdot x_i + b) \geq 1 \quad \forall i
\end{aligned}$
dove $ w $ è il vettore dei pesi, $ b $ è il bias, $ x_i $ è il vettore dei punti e $ y_i $ indica la classe.
Successivamente al calcolo del lagrangiano, si ottiene il seguente problema :
\begin{cases}
\min_{\boldsymbol{\alpha}} \frac{1}{2}\boldsymbol{\alpha}^T Q \boldsymbol{\alpha} - \sum_{i=1}^n \alpha_i\\
\sum_{i=1}^n \alpha_i y_i = 0\\
\alpha_i \geq 0\,,\quad \forall \ i=1,\ldots ,n
\end{cases},,
dove $\boldsymbol{\alpha}\in\mathbb{R}^n$ e la matrice $Q\in\mathbb{R}^{n\times n}$ è tale che
$Q = \left(q_{i,j}\right)_{i,j=1,\ldots ,n} = \left( y_iy_j\boldsymbol{x}_i^T\boldsymbol{x}_j \right)_{i,j=1,\ldots ,n}\,$.
A questo punto possiamo scrivere la funzione di decisione, ovvero quella che ci indica la classe di ciascun nuovo vettore, in questo modo:
${f}(x)=\text{sgn}(w · x + b )=\text{sgn}(\sum_{i=1}^n a_i y_i(x_i \cdot x)+b)$
con $b=y_i - \sum_{j=1}^n a_j y_j (x_j \cdot x_i)$
Caso non separabile
Nel caso non separabile, non esiste un iperpiano che possa separare perfettamente i punti di training.
In questo caso, le SVM cercano di trovare l'iperpiano che massimizza il margine e allo stesso tempo minimizza l'errore di classificazione sui punti di training.
Per ottenere ciò, le SVM introducono le variabili slack, che misurano la quantità di violazione del margine consentita per ciascun punto dato.
Il problema Primale diventa quindi
\begin{cases} \min_{\boldsymbol{w}} \frac{1}{2}\boldsymbol{w}^T\boldsymbol{w} + C \sum_{i=1}^n \xi_i\\ y_i (\boldsymbol{w}^T\boldsymbol{x}_i + b) \geq 1 - \xi_i \,,\quad & \forall \ i=1,\ldots ,n\\ \xi_i\geq 0\,,\quad & \forall \ i=1,\ldots ,n\\ \end{cases},,
Il parametro $C\in\mathbb{R}^+$ è un parametro di regolarizzazione che caratterizza il rilassamento delle condizioni per il margine:
- $C\rightarrow 0$ aumenta la "morbidezza" del margine, permettendo ai vettori $\boldsymbol{x}_i$ di superarlo illimitatamente;
- $C\rightarrow +\infty$ aumenta la "durezza" del margine, permettendo ai vettori $\boldsymbol{x}_i$ di superarlo impercettibilmente;
#calcolo la precisione del metodo SVM hard
C_hard = 1e10
loss = 'squared_hinge'
dual = False
lsvm_hard = LinearSVC(C=C_hard, loss=loss, dual=dual, random_state=random_seed)
lsvm_hard_pca = LinearSVC(C=C_hard, loss=loss, dual=dual, random_state=random_seed)
lsvm_hard.fit(X_train, y_train)
lsvm_hard_pca.fit(pca_adult, y_train)
df_lsvm_hard = pd.DataFrame({'acc. (no PCA)': [np.round(lsvm_hard.score(X_train, y_train)*100,decimals=2), np.round(lsvm_hard.score(X_test, y_test)*100, decimals=2)],
'acc. (PCA)': [np.round(lsvm_hard_pca.score(pca_adult, y_train)*100,decimals=2), np.round(lsvm_hard_pca.score(pca.transform(scaler_X.transform(X_test)), y_test)*100, decimals=2)]},
index=['training', 'test'])
display(df_lsvm_hard)
| acc. (no PCA) | acc. (PCA) | |
|---|---|---|
| training | 79.44 | 83.79 |
| test | 78.93 | 84.09 |
#calcolo la precisione del metodo SVM soft
C_soft = 1e1
loss = 'squared_hinge'
dual = False
lsvm_soft = LinearSVC(C=C_soft, loss=loss, dual=dual, random_state=random_seed)
lsvm_soft_pca = LinearSVC(C=C_soft, loss=loss, dual=dual, random_state=random_seed)
lsvm_soft.fit(X_train, y_train)
lsvm_soft_pca.fit(pca_adult, y_train)
df_lsvm_soft = pd.DataFrame({'acc. (no PCA)': [np.round(lsvm_soft.score(X_train, y_train)*100,decimals=2), np.round(lsvm_soft.score(X_test, y_test)*100, decimals=2)],
'acc. (PCA)': [np.round(lsvm_soft_pca.score(pca_adult, y_train)*100,decimals=2), np.round(lsvm_soft_pca.score(pca.transform(scaler_X.transform(X_test)), y_test)*100, decimals=2)]},
index=['training', 'test'])
display(df_lsvm_soft)
| acc. (no PCA) | acc. (PCA) | |
|---|---|---|
| training | 79.44 | 83.79 |
| test | 78.93 | 84.09 |
Notiamo che sia col metodo "hard" sia col metodo "soft" otteniamo la stessa accuratezza.
Mentre invece se applichiamo il metodo SVM al training sample già ridotto dalla PCA l'accuratezza aumenta circa del 5%
SVM Non lineare¶
Introduzione
I Support Vector Machines (SVM) non lineari sono un'estensione degli SVM lineari che possono separare dati non linearmente separabili utilizzando il "Kernel's trick".
Questo metodo trasforma i dati in uno spazio di dimensioni superiori dove un iperpiano lineare può separare le classi.
Procedimento
Sia $X$ lo spazio dei dati, definiamo un funzionale $ \phi : X \rightarrow {H} $ dove ${H}$ è uno spazio di Hilbert dove ${H}$.
Definiamo ora un Kernel come $K(x , x')=\langle\phi(x),\phi(x')\rangle$ dove $x,x'$ sono vettori.
Mappatura in uno Spazio di Dimensioni Superiori:
$$\phi :\mathbb{R}^n\rightarrow \mathbb{R}^m\,, \quad \text{con }m>n\,,$$ I dati originali vengono mappati in uno spazio di Hilbert di dimensioni superiori utilizzando il funzionale $\phi(x) $.
In questo spazio, è più probabile che i dati siano linearmente separabili.
Siano ${X}=\{\boldsymbol{x}_1,\ldots ,\boldsymbol{x}_n\}$ e ${Y}=\{\boldsymbol{y}_1,\ldots ,\boldsymbol{y}_n\}$ i vettori dei dati e delle classi, allora il vettore $X$ diventa $X^\phi=\{ \phi(x_1), \ldots , \phi(x_n)\}=\{ \boldsymbol{\varphi}_1,\ldots ,\boldsymbol{\varphi}_n\}$.Problemi con la soluzione precedente
Come già mostrato per SVM lineare bisogna risolvere lo stesso problema di ottimizzazione ma sta volta utilizzando $X^\phi$, le soluzioni sono:
$$\boldsymbol{w}^* = \sum_{i=1}^n y_i\alpha_i^* \boldsymbol{\varphi}_i$$
$$b^* = \frac{1}{n} \sum_{i=1}^n (y_i - \langle\boldsymbol{w}^*, \boldsymbol{\varphi}_i\rangle)\,$$
Ora possiamo scrivere la nostra funzione di decisione:
\begin{aligned} \ {f}(x) &= \mathrm{sign}\left( \langle \boldsymbol{w}^* , \phi(\boldsymbol{x})\rangle + b^*\right) = \\ & = \mathrm{sign}\left( \sum_{i=1}^n y_i\alpha_i^* \langle \boldsymbol{\varphi}_i , \boldsymbol{\varphi}\rangle + \frac{1}{n} \sum_{i=1}^n \left(y_i - \sum_{j=1}^n y_j\alpha_j^* \langle\boldsymbol{\varphi}_j , \boldsymbol{\varphi}_i\rangle \right) \right)\, \end{aligned}
Nella formula sopra, il grosso dei calcoli è richiesto per i prodotti scalari $\langle \boldsymbol{\varphi}_i , \boldsymbol{\varphi}\rangle$ e $\langle\boldsymbol{\varphi}_j ,\boldsymbol{\varphi}_i\rangle$, che vanno svolti in $\mathbb{R}^m$, dove presumibilmente $m\gg n$ (teoricamente anche $m=\infty$).
Per questo dobbiamo trovare un modo per semplificare i calcoli.Scelta del kernel
Come definito in precedenza, un kernel è una funzione che calcola il prodotto scalare tra due vettori nello spazio di dimensioni superiori senza calcolare esplicitamente $\phi(x)$.
Mostriamo ora alcuni esempi di kernel più utilizzati(noi faremo un esempio utilizzandone due distinti):Kernel Lineare:
$K(x_i, x_j) = \langle x_i , x_j \rangle$
Kernel Polinomiale:
$K(x_i, x_j) = (\gamma\langle x_i , x_j \rangle+ c)^d$
Kernel Radiale di Base (RBF):
$K(x_i, x_j) = \exp\left(-\gamma \langle x_i , x_j \rangle^2\right)$
Kernel Sigmoide:
$K(x_i, x_j) = \tanh(\gamma \langle x_i , x_j \rangle + c)$
Una volta scelto il kernell più adatto lo si sostituisce al prodotto scalare e si ottiene la funzione di decisione semplificata.
#testo il primo algoritmo con un kernel sigmoid, tempo medio 29 secondi
ker_tanh = 'sigmoid'
gamma_tanh = 'auto'
C = 100
svm_tanh_1 = SVC(C=C, kernel=ker_tanh, gamma=gamma_tanh, random_state=random_seed)
svm_tanh_2 = SVC(C=C, kernel=ker_tanh, gamma=gamma_tanh, random_state=random_seed)
svm_tanh_1.fit(X_train, y_train)
svm_tanh_2.fit(pca_adult, y_train)
df_svm_tanh = pd.DataFrame({'acc. (NO PCA)': [np.round(svm_tanh_1.score(X_train, y_train)*100, decimals=2), np.round(svm_tanh_1.score(X_test, y_test)*100,decimals=2)],
'acc. (PCA)': [np.round(svm_tanh_2.score(pca_adult, y_train)*100,decimals=2), np.round(svm_tanh_2.score(pca.transform(scaler_X.transform(X_test)), y_test)*100,decimals=2)]},
index=['training', 'test'])
df_styled=df_svm_tanh.style.set_caption('SVM SIGMOID')
display(df_styled)
| acc. (NO PCA) | acc. (PCA) | |
|---|---|---|
| training | 75.370000 | 77.210000 |
| test | 74.840000 | 77.600000 |
#testo il secondo algoritmo con un kernel rbf, tempo medio 3 minuti 38 secondi
ker_rbf = 'rbf'
gamma_rbf = 'auto'
C = 100
svm_rbf_1 = SVC(C=C, kernel=ker_rbf, gamma=gamma_rbf, random_state=random_seed)
svm_rbf_2 = SVC(C=C, kernel=ker_rbf, gamma=gamma_rbf, random_state=random_seed)
svm_rbf_1.fit(X_train, y_train)
svm_rbf_2.fit(pca_adult, y_train)
df_svm_rbf = pd.DataFrame({'acc. (NO PCA)': [np.round(svm_rbf_1.score(X_train, y_train)*100, decimals=2), np.round(svm_rbf_1.score(X_test, y_test)*100,decimals=2)],
'acc. (PCA)': [np.round(svm_rbf_2.score(pca_adult, y_train)*100,decimals=2), np.round(svm_rbf_2.score(pca.transform(scaler_X.transform(X_test)), y_test)*100,decimals=2)]},
index=['training', 'test'])
df_styled=df_svm_rbf.style.set_caption('SVM RBF')
display(df_styled)
| acc. (NO PCA) | acc. (PCA) | |
|---|---|---|
| training | 99.990000 | 90.110000 |
| test | 74.190000 | 82.740000 |
Per la scelta del kernel migliore bisogna cercare di fare un bilanciamento tra precisione ottenuta e tempo di calcolo.
#ora disegno i grafici per il kernel rbf, considerando solo le prime due componenti principali tempo medio 2 minuti 8 secondi
ker_rbf = 'rbf'
gamma_rbf = 'auto'
C = 100
svm = SVC(C=C, kernel=ker_rbf, gamma=gamma_rbf, random_state=random_seed)
svm.fit(pca_adult[:, :2], y_train)
df_svm_rbf = pd.DataFrame({'acc. (PCA)': [np.round(svm.score(pca_adult[:, :2], y_train)*100,decimals=2), np.round(svm.score(pca.transform(scaler_X.transform(X_test))[:, :2], y_test)*100,decimals=2)]},
index=['training', 'test'])
display(df_svm_rbf)
minx=min(pca_adult[:, 0].min()-1, pca_adult[:, 1].min()-1)
maxx=max(pca_adult[:, 0].max()+1, pca_adult[:, 1].max()+1)
mesh_points = 200
xx1, xx2 = np.meshgrid(np.linspace(minx, maxx, mesh_points), np.linspace(minx, maxx, mesh_points))
XX = np.concatenate([xx1.reshape(mesh_points ** 2, 1), xx2.reshape(mesh_points ** 2, 1)], axis=1)
ZZ_1 = svm.decision_function(XX)
zz_1 = ZZ_1.reshape(mesh_points, mesh_points)
plt.figure(figsize=(12, 12))
plt.scatter(pca_adult[:, 0], pca_adult[:, 1], c=y_train,marker='x', cmap=cmap, alpha=0.25)
plt.scatter(pca.transform(scaler_X.transform(X_test))[:, 0], pca.transform(scaler_X.transform(X_test))[:, 1], c=y_test,marker='o', cmap=cmap, alpha=0.25)
plt.contour(xx1, xx2, zz_1, [0], colors='red')
plt.contour(xx1,xx2,zz_1,[-1,1],colors=['green','yellow'])
plt.title('SVM RBF - DECISION BOUNDARY')
| acc. (PCA) | |
|---|---|
| training | 81.97 |
| test | 81.77 |
Text(0.5, 1.0, 'SVM RBF - DECISION BOUNDARY')
y_pred = svm_rbf_1.predict(X_test)
y_pred_pca = svm_rbf_2.predict(pca.transform(scaler_X.transform(X_test)))
precision = np.round(precision_score(y_test, y_pred)*100,decimals=2)
recall = np.round(recall_score(y_test, y_pred)*100,decimals=2)
fbeta = np.round(fbeta_score(y_test, y_pred, beta=1.0)*100,decimals=2)
precision_pca = np.round(precision_score(y_test, y_pred_pca)*100,decimals=2)
recall_pca = np.round(recall_score(y_test, y_pred_pca)*100,decimals=2)
fbeta_pca = np.round(fbeta_score(y_test, y_pred_pca, beta=1.0)*100,decimals=2)
df_scores = pd.DataFrame({'Precision': [precision, precision_pca], 'Recall': [recall, recall_pca], 'F1': [fbeta, fbeta_pca]}, index=['NO PCA', 'PCA'])
df_styled = df_scores.style.set_caption('SVM RBF - SCORES')
display(df_styled)
| Precision | Recall | F1 | |
|---|---|---|---|
| NO PCA | 43.360000 | 8.430000 | 14.120000 |
| PCA | 67.360000 | 60.910000 | 63.970000 |